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ABSTRACT 

High-energy particles stream during coronal mass ejections or flares through the plasma of the solar wind. This causes instabilities, 
which lead to wave growth at specific resonant wave numbers, especially within shock regions. These amplified wave modes influence 
the turbulent scattering process significantly. In this paper, results of particle transport and scattering in turbulent plasmas with excited 
wave modes are presented. The method used is a hybrid simulation code, which treats the heliospheric turbulence by an incompressible 
magnetohydrodynamic approach separately from a kinetic particle description. Furthermore, a semi-analytical model using quasilinear 
theory (QUT) is compared to the numerical results. This paper aims at a more fundamental understanding and interpretation of the 
pitch-angle scattering coefficients. 

Our calculations show a good agreement of particle simulations and the QUT for broad-band turbulent spectra; for higher turbulence 
levels and particle beam driven plasmas, the QUT approximation gets worse. Especially the resonance gap at yu = poses a well- 
known problem for QLT for steep turbulence spectra, whereas test-particle computations show no problems for the particles to scatter 
across this region. The reason is that the sharp resonant wave-particle interactions in QUT are an oversimplification of the broader 
resonances in test-particle calculations, which result from nonlinear effects not included in the QUT. We emphasise the importance of 
these results for both numerical simulations and analytical particle transport approaches, especially the validity of the QUT. 

Key words. Magnetohydrodynamics (MHD) - Turbulence - Sun: coronal mass ejections (CMEs) - Sun: particle emission - acceler- 
ation of particles 



1. Introduction 

Coronal mass ejections (CME) are believed to be an important 
source of solar energetic particles (SEP). The acceleration mech- 
anism taking place is the diffusive shock acceleration process 
([Reames 1999 ) One cornerstone in the theory behind this mech- 
anism is understanding particle scattering of SEP in the back- 
ground plasma. The general approach used for this purpose is 
the quasilinear theory (QLT) ( Jokipii|1966 1 

The CME may be represented as a system of a thermal 
plasma with a turbulent energy spectrum and a nonthermal pro- 
ton and electron component. These constituents are interact- 
ing with each other in a nonlinear way. In order to understand 
this system, it is necessary to model wave excitation by parti- 
cle beams, turbulent energy transfer, and the transport of high- 
energy particles. Models have been developed by, e.g. Vainio & 
|Laitinen| ( [2007l l and |Ng & Reames| ( [2008] l While these models 
are capable of describing the interaction of particles and waves, 
the details of the specific processes are implemented in a sim- 
ple way, basically relying on QLT or its simple extensions as the 
correct description of particle scattering and ignoring oblique or 
perpendicular wave modes completely. Lange & Spanier ( 2012| l 
did numerical simulations of the turbulent cascading process of 
energy injected via proton beams. The basic concept of this pa- 
per is to take first steps towards a consistent description of parti- 
cle transport within these plasma configurations and investigate 
how well the QLT describes scattering in such turbulence. 

The general problem in the derivation of transport coeffi- 
cients is a gap between theoretical and numerical approaches. 
Numerical simulations have used artificial turbulence spectra 



( Qin et al.|2002| l! while analytical calculations use any possible 
spectrum. However, the applicability of these calculations is un- 
clear. This work combines the previous turbulence simulations 
by Lange & Spanier (2012 ) with a particle tracking algorithm 
derived from S panier & Wisniewski (201 1 ) The results are com- 
pared with QLT calculations using an identical input spectrum. 

This approach allows a comparison of QLT and numerical 
simulations. The latter are assumed to be valid for a wide range 
of turbulent spectra and particle energies. Therefore, a determi- 
nation of the validity of QLT is possible. Novel visualization 
and comparison methods are used to probe nonlinear and time- 
dependent scattering processes. Thus, a more fundamental un- 
derstanding of the scattering process is achieved. 



2. Theory 

2.1. Numerical Model 

The region of the heliosphere we are interested in is within the 
weak turbulence regime, with magnetic field fluctuations defined 

as 



6B = B-B, 



0, 



(1) 



where {6B} = is the mean value of the fluctuations, which 
leads to {B) = Bq. In that sense, weak turbulence is connected 
to small fluctuation amplitudes compared to the magnetic back- 
ground field Bq. In terms of wave turbulence, this means that 
the linear Alfven wave timescale is shorter than the intrinsically 
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nonlinear time. It is observed that the solar wind magnetic fluc- 
tuations decrease as 5B^ oc r"^, while the background field de- 



creases as B}. oc r (Bavassano et al. 1982; Bruno & Carbone 
|2005)l| Consequently the dJJ/JJo ratio withm the heliosphere is 
increasing with distance from the Sun ( Hollweg et al.|2010) l 

The magnetic background field Bq is defined towards the 
z-direction within our simulations and hence is also denoted 
as The parallel direction is, therefore, defined as the z- 

direction and the x- and y-direction are the perpendicular direc- 
tions. For symmetry reasons, there will be no further distinction 
between the two perpendicular directions, and all plots show val- 
ues averaged over the azimuthal angle in cylindrical coordinates 
of the x-y plane. 

High Reynolds numbers in combination with massive energy 
injection, as seen in, e.g. the solar wind, are strong indicators of 
a highly turbulent state. This means that most parts of the helio- 
sphere are dominated by turbulent evolution and consequently 
energy cascading to smaller spatial scales. In situ measurements 
of the energy spectrum ( Tu et al.]|198 9 ) are in agreement with 
this fact. 

To simulate conditions within the turbulent heliospheric 
plasma with particle transport and scattering, the research group 
at the University of Wiirzburg has developed a hybrid simulation 
code, GisMO. It is divided into two parts, the magnetohydrody- 
namic (MHD) algorithm Gismo-MHD and the kinetic particle 
transport simulation tool Gismo-Particles. 

Gismo-MHD is an incompressible pseudospectral MHD 
code that is fully parallelised and capable of efficiently using 
massive computing clusters. For a more detailed description we 
refer to Lange & Spanier (2012 1 A brief overview is presented 



below. The basis of the simulation software is to solve the fol- 
lowing set of incompressible MHD equations: 



du 

Ik 
db 

¥ 

V ■ u 



bVb-uWu-WP + v,V^''u 



bVu-uVb + Tj^'-'^b 


0, 



(2) 



where b - B/ yjA-ng is the normalised magnetic field, u is the 
fluid velocity, and g is the constant mass density. The dissipative 
terms due to viscous and Ohmic dissipation are denoted by v,, 
and rj. In the following we introduce the parameter v as a global 
diffusivity with rj - = v. The approach is valid for mag- 
netic Prandtl number of order unity, which is applicable within 
the regime of Alfven wave turbulence where an equipartition be- 
tween magnetic and kinetic energy can be assumed ( Bershadskri 
|2002[ |Bigot et "aL]|2008[ ) The artificial increase of the dissipa- 
tion by the hyperdiffusivity parameter h is often used in spectral 
methods. The pressure term Vf fulfils the closure condition for 
incompressibility ( |Maron & Goldreich|2001[)j 



V^f = Vft : Vft - Vh : Vh. 



(3) 



In the incompressible regime of a magnetised plasma, the 
MHD turbulence consists of only two types of waves, which 
transport energy along the parallel direction: the so-called 
pseudo and shear Alfven waves. The former are the incompress- 
ible limit of slow magnetosonic waves and play a minor role 
within anisotropic turbulence ( Maron & Goldrei"ch||2001| l The 
pseudo Alfven waves polarisation vector is in the plane spanned 
by the wave vector k and Bq. The shear waves are transver- 
sal modes with a polarisation vector perpendicular to the k - 



Bq plane. They are circularly polarised for parallel propagat- 
ing waves. Both species exhibit the dispersion relation oj^ - 
{vAk\\f'. We note that the shear mode seems to be dominant be- 
cause pseudo waves are heavily damped by the Barnes damping 
process within weakly turbulent regimes ([Sridhar & Goldreich 



1994) The damping weakens in strong turbulence, but according 
to Goldr eich & Sridhar| ( |1995] l the wave generation of pseudo 
Alfvenic wave modes is only possible via three-wave interac- 
tions by two shear wave modes. Barnes damping is important for 
high-j6 plasmas. Since this is not fulfilled for the solar corona, the 
role of pseudo Alfven waves should not be ignored. However, 
the solar wind and especially its fast component is strongly dom- 



inated by Alfven waves (Bruno & Carbone 2005 j Howes et al. 



|2012 liand therefore incompressible MHD is well applicable. 

Because of this model, a suitable description of the system is 
achieved by the use of Alfvenic waves moving either forwards 
or backwards. Therefore the Elsasser variables (El sasser||1950|l| 
are introduced: 



w -V + b - vaBw 
-V - b + v^ey. 



(4) 



Applying this definition to Eqs. (|2]) and transformation to the 
Fourier space yields 



{d, - va^z) w„ = 



i kak/jky / ^ 



+ w„ WZ 



(5) 



where the tilde-notation represents quantities in the Fourier 
space. This is the final set of equations which is iterated by 
Gismo-MHD. 

Obviously, the nonlinearities of Eqs. (|5]l that describe the tur- 
bulent behaviour of the MHD plasma cannot be solved in Fourier 
space in a straightforward fashion. Hence the main numerical 
load is the transfomation between real and wavenumber space 
for each iterative step. For this purpose we used the P3DFFT al- 
gorithm, which is an MPI-parallelised fast Fourier transfomation 
based on FFTW3 ( Pekurovsky 2011 ) 

A basic problem of spectral methods that use discrete Fourier 
transformation is the aliasing effect. Because of discrete sam- 
pling in the wavenumber space, high fc-values exhibit errors that 
depend on the structure of the real space fields. Therefore we 
used zero padding, which is also referred to as Orszag's 2/3 rule, 
i.e. 2/3 of the wavenumbers below the Nyquist frequency have to 
be truncated to achieve maximum anti-aliasing, hence reducing 
the Fourier space resolution to 1/3 of the original wavenumber 
range (Orszag 1971 ) This process is repeated for each step, im- 
mediately before calculating the nonlinearities and, accordingly, 
calculating the right-hand side (RHS) of the MHD equations. 
Consequently, the change in the antialiasing-range of one MHD 
step is physically correct, but not the long-term evolution. 

The code Gismo is capable of using different foward-in-time 
schemes, namely, Euler and Runge-Kutta (RK) second as well 
as fourth order. All the simulations in this paper have been pe- 
formed using RK-4. 

The Alfven wave generation mechanism by SEPs is not in- 
vestigated in detail here. The driving mechanism assumed for 
our simulations is the streaming instability. The estimation of the 
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wave growth rate is described in |Vainio| ( |2003| l The streaming 
instabihty is caused by energetic proton scattering off the coro- 
nal Alfven waves. During the the scattering process the particle 
changes its pitch-angle cosine by A//, while its wave-frame en- 
ergy remains constant. Thus the particles' energy in the plasma 
frame is changed by VApAfi, where p is the wave-frame par- 
ticle momentum. Accordingly, the Alfven wave energy in the 
plasma frame changes due to energy conservation. Another im- 
portant instability in the solar corona is the electrostatic insta- 
bility, which is caused by an electron current and streaming 
ions. Ion acoustic waves would be generated by this process. 
However, for the growth rate of these modes a sufficiently high 
ratio Tg/Ti » 1 of the electron and ion temperatures is crucial. 
Observations and simulations in the vicinity of three solar radii 



indicate temperature ratios of the order of unity (Landi 2007 



Jin et al. 2012) In this regime the ion acoustic waves are also 
efficiently suppressed by Landau damping. For further reading 



about the streaming instability we refer to [Gary ( 1993 1 

The numerical tool for the kinetic simulation of charged test 
particles is Gismo-Particles. Like Gismo-MHD it is fully paral- 
lelised and capable of using massive computing clusters. Its core 
is the calculation of the relativistic Lorentz force 

^ jv ^ —[cE(x,t) + vxB(x,t)], (6) 
df mc 

which acts on the particles at position x through the electric E 
and magnetic fields B of the plasma. The particle velocities are 
denoted by v, c is the speed of light and y represents the Lorentz 
factor The Gismo-Particles tool is capable of simulating test 
particles using physical masses and charges of electrons and 
ions. The test particles do not influence the background plasma, 
i.e. self-consistent backreactions to the plasma are neglected. A 
charged particle within a magnetic field will perform a gyromo- 
tion with the frequency 



Q 



ZeB 
ymc 



and the Larmor radius 



(7) 



(8) 



A useful numerical approach for solving Eq. (|6]l for gyrating par- 
ticles is the implicit scheme of the Borispush. The basic idea was 
given by Boris ( 1970[ l where the iterations of the Lorentz force 
are separated in two partial steps. First, the particles are accel- 
erated by the electric field within a half time step. Second, the 
gyromotion of the particles, which is caused by the magnetic 
field, is calculated. After that, the electric fields acts again for 
another half time step to complete the iteration. This approach 
leads to a discretisation of the Lorentz force with the following 
set of equations: 



lowing steps were used: 

v' 
hi 
v + 



q B 



-At" 



2 y" m c 

2hi 
\ + h\ hi 
V " -H v' X hi. 



(12) 
(13) 
(14) 
(15) 



where the auxiliary vectors h were introduced to calculate v ^ 
(Birdsall & Langdon 2005 ) By using these equations, the new 
velocity and location for each particle are then 



X + 



qE At 
2 m 

y(y r+Ar/2-) 



Af" 



(16) 



(17) 



The advantage of the Borispush is the very high numerical sta- 
bility. The particles are assumed to undergo gyromotions, hence 
the particle orbits themselves are stable for an arbitrary time dis- 
cretisation. Even in the limit of Af""'" » Q"', the particle orbit is 
stable but converges to an adiabatic drift motion. The limitation 
of this method is the correct resolution of the Larmor radius r^. 
If the timestep chosen is too large, this leads to a big deviation 
from the analytical r^. 

To provide a correct resolution of the gyromotion, Gismo- 
Particles performs a couple of single-particle gyrations with 
the used parameters during each start-up phase. The Larmor ra- 
dius is then measured using a three point-circle approximation. 
If the deviation to the analytical value is above a given accu- 
racy threshold, the gyrosimulation is performed with a smaller 
timestep Af""™. This procedure repeats until the accuracy condi- 
tion is fulfilled. Specifically, in our simulations an accuracy of 
the order of |r£ - rmeasui-edl/'"L ~ 10"^ was used. This ensures that 
in each simulation the time discretisation is chosen correctly. Of 
course, the discretisation of spatial grid must resolve as well. 
This means that twice the Larmor radius must be much larger 
than a grid cell and smaller than the whole system length. 

Ultrarelativistic particle speeds present limitation to the 
method of the Borispush. In this case the conservation of en- 
ergy would be violated, since the ideal ohmic law is not fulfilled 
anymore. Beyond Lorentz factors of y x 10^, fictitious forces 
start to act and this approach is no longer applicable ( Vay 2008]! 



Both parts of Gismo are calculated for each step. After 
iterating the Elsasser MHD fields w'^, they are transformed 
into the physical, electric, and magnetic fields which are trans- 
ferred to Gismo-Particles. Then the Borispush is performed. 
Each particle responds to its local fields, which are calculated 
by an averaging method via three-dimensional splines (Spanier] 
& Wisniewski 2011 ; Wisniewski et al. 2012) Periodic spatial 
boundary conditions were used; thus the number of particles re- 
mained constant in each simulation. 



1 



2 ym 



q E At"""" 
2 m 

(v^ +v')xB 



At 



, t+&t/2 _..+ , qE At" 



2 m 



(9) 
(10) 

(11) 



where quantities with Af""™/2 denote the discrete half-time 
steps. To solve this set of equations with respect to v ^, the fol- 



2.2. Statistical description 

The general statistical approach of particle transport by the rela- 
tivistic Vlasov equation, 

^ + A . ^ + i (cE(x, + V X B(x, 0) ■ ^ = St(x,p, t), 
at my ox c op 

(18) 

describes the development of a particle distribution fr of species 
T under the influence of the Lorentz force (|Schlickeiser|2002^ 
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The momentum is denoted by p and the term on the RHS 5 7- is 
a source term. It is useful to transform this equation into guiding 
centre coordinates, since we are not interested in the actual posi- 
tion of the gyration, but in the averaged position of the particle. 
The Vlasov equation then reads 

dpT OFt OFt Id t 

+ - + ^ ^{{p'gx^dfT)) = SriX^, f). 



dt 



dZ 



(19) 



where Fj - (fr) is the expectation value of fr and X^r are the 
guiding centre coordinates with the angle in the perpendicular 
plane and the pitch-angle cosine yu. These coordinates are given 
by 



= p -^^r^^cos 

Py = P 



fi^ sin ^ 



Pz= P^i 



x = X- —pJl -H^ sine 
qBo y 



(20) 



The generalised forces gx^ represent the effects of the electro- 
magnetic fields and are basically the time derivatives of the de- 
noted variables, e.g. gx^ - Xa-. This equation is in general ana- 
lytically unsolvable. 

A common approach to describing particle transport analyt- 
ically is the QLT, which was first suggested by Jokipii ( |1966| l 
in the context of energetic charged particle transport in turbulent 
magnetic fields. Its core is the assumption of unperturbed parti- 
cle orbits. This implies that the fluctuation amplitudes are small, 
leading to a quasilinear system. The Vlasov equation Eq. ( [T9] l 
would then simplify to 

OFt OFt OFt 



1 d 

P^dXa 



I OFt 
dXr 



Jo 



^gxiX^,s)) 



(21) 

where the method of characteristics was applied and the hat sym- 
bol represents quantities along the characteristics ( Schlickeiser] 
|2002| ) This equation is known as the Fokker-Planck equation 
with the Fokker-Planck coefficients D 



One of the most in- 
teresting variables is the pitch-angle diffusion coefficient D^^. It 
describes the pitch-angle scattering of the particle and is conse- 
quently connected to the scattering mean free path, which can be 
evaluated by the observable angle distribution and Monte Carlo 
simulations ( Agueda et al. 2009 ) Here, scattering means a reso- 
nant wave-particle interaction of nth order which fulfils the con- 
dition 



k\\ V\\ — CO + n D. — 0, n € 



(22) 



( Schlickeiser|| 1 989 1' where cj is the wave frequency and fcy its 
parallel wavenumber. Q is again the particle gyrofrequency and 
V|| its parallel velocity component. Whether a particle with vy 
interacts resonantly with a wave with Ary is determined by the 
polarisation ( jKennel & Engelmann| [l966) Because our MHD 
model consists of pseudo and shear Alfven waves only, certain 



values of n can be connected to the different types of waves. The 
Cherenkov resonance, n=0, is generated only by waves with a 
compressive magnetic field {6R ■ Bq 0), i.e. pseudo Alfven 
waves. Purely parallel waves can contribute only to the n = ±1 
resonance, and resonances with \n\ > 1 occur only for waves 
with nonvanishing perpendicular wavenumber, i.e. with oblique 
Alfven waves. For a mathematical treatment of the resonant in- 
teractions, see Appendix [B| 

For further reading and a detailed insight of the derivation 
of Eq. 



(21 



we refer to Schlickeiser ( 2002]) A serious problem 
of the QLT is the limited applicability. The approximation of 
small fluctuations only holds for weak turbulent systems, where 
6B/B() <K 1. This is important for the local field which acts on the 
invidual particles. For instance, a strong turbulent region would 
change effectively the direction of the mean magnetic field and 
consequently the gyromotion of the particles. Hence the assump- 
tion of unperturbed orbits would be invalid. Another problem of 
the QLT is the inadequate description of particles propagating 
perpendicular to the mean magnetic field, fi x Q. However, re- 
garding Eq. ( 22 1 a very small parallel particle velocity will gen- 
erate a singularity for vy — > 0. One aspect of this paper will 
concentrate on the applicability of the QLT to describe solar en- 
ergetic particle transport. 

2.3. The Fokker-Planck coefficient D^jf^ 

We will now focus on the pitch-angle scattering coefficient in 
more detail. The basic idea is to compare the simulation results 
of GisMO with the analytical approach. 

Our simulations provide us with all the information of the 
trajectory and velocity of each individual particle within a turbu- 
lent plasma. Particle scattering will be considered here in terms 
of the pitch-angle cosine 



fi - cos(q') 



V B 

W\\B\' 



(23) 



which describes the orientation of the particles' velocity vector 
with respect to the local magnetic field. In the approximation of 
the QLT, there are different ways to calculate by using the 
absolute change of yU during a certain time interval for a single 
particle. We used the definition 



lim 



(AmY t»'o (Apr 



2 At 



2 At 



(24) 



where the time interval At - t - to compared to an initial state 
at to is assumed to be large, i.e. the time evolution t has to be 
sufficient to develop resonant interactions. The necessary time 
development is discussed in the result section. Instead of using 
the change of the pitch-angle cosine A^ - ji^- po, the definition 
with the change of the pitch-angle Aa - Ue - ao itself leads to 



(Aaf 
2 At 



D 



1 



/2' 



(25) 



which represents the scattering frequency. 

For the analytical approach we use the magnetostatic limit 
for the calculation of D^^ for pseudo 



2ff(i -//) 



d^k7T6{k\\v\\ + nQ) X 



n#0 

2 



(26) 
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and shear Alfven waves 

2n\l - "2 



«^«(v'±fc±/n) 



P„,A(k), (27) 



where f „ is the jcx-component of the magnetic field fluctuation 
tensor, which represents the turbulence power spectrum. These 
equations are based on Schlickeiser (2002) and a more detailed 
derivation is presented in AppendixjB] We note that a proper de- 
scription of the Cherenkov resonance in Eq. 26 is not possible 
in this formulation. Using n = yields an undefined mathemati- 
cal expression. These coefficients were semi-analytically solved 
with the set of equations given in Appendix [C] We note that the 
approach of the magnetostatic QLT (SQLT) assumes the wave 
frequency to be small compared to the gyro frequency, w <s Q. 

3. Simulation setup 

We focus on the environmental conditions in the solar corona at 
a distance of three solar radii. The magnetic background field 
in this region is approximately Bq - 0.174 G, which yields, 
assuming a particle d ensity of 10^ cm~ ^, an Alfven speed of 
va - 1-2 ■ 10** cm s"' ( Vainio et al.|2003 i This region is of spe- 
cial interest because particle acceleration by CME-driven shocks 
is strongest there. 

The outer length scale of the simulated system is Lscaie - 
3.4 • lO^cm. This value is estimated using the growth rate from 
|Vainior|2003l ) 



nyQ. 
2n„VA 



A^pv^i\k\6m--)f, (28) 

V 



with proton cyclotron frequency Q., proton speed v, the Lorentz 
factor y, the proton number density Hp, as well as ji, the pitch- 
angle cosine, and / the proton distribution function. As the reso- 
nance condition, Eq. ( 22 1 is applied, assuming that the particles 
are streaming along the field. 

A remark on the notation: the wavenumber is in general de- 
fined as A; = {2nj)IL, where j stands for the numerical grid 
position. For simplicity we used the normalised wavenumbers 
k' - kL- {2n ■ j) throughout. 

We generated peaks at k'^^ -2n-% and k'^^ - 2n ■ 24, which can 
be excited by streaming protons with energies of £ a; 64 MeV 
and E ^ 7 MeV respectively. Using the resonance condition, this 
leads to a length scale of 



'scale 



27r /' , 

ynipcv a 10 cm. 

cBq 



(29) 



On this scale the magnetic background field is assumed to be 
homogeneous. The spatial resolution is 256"' gridpoints, result- 
ing in 128"' points in k-space of which \k'\ - 2ji[Q ■ ■ ■ 42] wave 
modes are active modes that remain unaffected by (anti)aliasing. 
The resistivity parameter was set to v = 1 in numerical units and 
the hyperdiffusivity coefficient ioh -2. 

The turbulence was simulated using an anisotropic driv- 
ing mechanism. Energy is, therefore, constantly injected into 
the simulation volume at the first 5 numerical wavenumbers 
in perpendicular {k'^ - 2ji[0 ■ • ■ 4]) and 15 in parallel direction 
(A;^ - 2n{Q ■ ■ ■ 14]). The reference frame is the solar wind frame. 

The anisotropy was chosen for two reasons: first, to mimic 
the preferential direction of the solar wind, where particles 



streaming radially away from the Sun form the Parker spiral. 
Consequently, these particles can deposit their energy in a paral- 
lel direction on different scales. This is mainly valid in the vicin- 
ity of the Sun, which we are interested in. Second, a slab compo- 
nent of solar wind turbulence is also observed at small scales in 
the parallel direction. To ensure turbulence evolution up to high 
parallel wavenumbers, the driving range was extended along the 
parallel axis. This is necessary because the parallel evolution is 
much weaker than the perpendicular one (jGoldreich & Sridhar 



1995 



1997 1 Even though this is primarily a technical aspect 
to ensure the extent of the spectrum to higher k\\, it is still in 
line with observations. An isotropic driver would not yield suffi- 
ciently turbulent modes at high k\\. 

The turbulence driving is performed by allocating an ampli- 
tude with a phase to the Elsasser fields within the Fourier space. 
The amplitude follows a power-law of \k['^^^ and is initialised 
using a normal distribution. The phase was randomly chosen 
between zero and 2:7r. These settings are divergence free and 
Hermitian symmetric. After this initialisation the values were 
scaled to the desired scenario, which in our case is a 6B/B() ra- 
tio of roughly 10"^. We note that both species, pseudo and shear 
Alfven waves, are excited by this type of turbulence driving, but 
as presented by |Maron & Goldreich ( 2001 1 the pseudo wave evo- 
lution is strongly suppressed. In this inertial range energy is in- 
jected every Q.03s (10 MHD steps), which leads to a saturated 
turbulence, an equilibrium between dissipation and injection. 

A Gaussian-distributed energy peak with purely parallel 
wavenumber k - ke\\ was injected within the saturated turbu- 
lence. For this purpose, two different wave modes were chosen. 
To investigate the physics of an SEP-event, a wavenumber of 
= 1.5 ■ 10"^ cm"' is used. This corresponds to a numerical 
wavenumber of 2:7r ■ 8, which is still within the driving range of 
the turbulence. To represent injection at smaller scales, a peak is 
set at = 4.4 ■ 10"^ cm"'. With k'^^ = 27r ■ 24, this wave mode is 
not in the driving range of the background turbulence. The injec- 
tion of energy in these modes was done gradually over a specific 
time interval. 

One important aspect of the calculations of D^p is the ap- 
plicability of the QLT. The ratio SB/Bq is, therefore, the limit- 
ing parameter (see Sect. 2.2 1. Consequently, it is of interest to 
explore peaks at either position under the influence of different 
magnetic background fields. For this purpose another turbulence 
simulation was performed which uses the same initial condi- 
tions, except for a higher mean field Bq = = 1-74 G and 
for stability reasons a higher dissipation rate v = 10cm^''s"'. 
In the following the smaller field simulation will be denoted by 
Bq = = 0.174 G. We note that the B„ case is an artificial 
scenario rather than representative of coronal conditions. Peaks 
in such high magnetic fields represent SEP energies up to 2.65 
GeV. In this case the SEPs cannot generate waves by streaming 
as the relativistic proton intensities are insufficient. 

To summarise, these four simulations with excited wave 
modes at k'^^ = 27r ■ 8 and k'^^ - 2-n ■ 24, both within turbulent 
fields governed by a strong and a weak Bq, are the starting point 
for the test particle simulations and the calculations of the scat- 
tering coefficient D^^. 

All of the test particle simulations by Gismo-Particles were 
performed with 10^ protons with an initial uniform distribution 
in IX and and random positions x. This amount has proven to 
be sufficient in test simulations for good statistics. These ini- 
tial conditions aim to provide a complete coverage of the phase 
space in yU for the test particles. Thus, we are not interested in the 
development of special distribution functions. 
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Table 1. Parameter setup for the particle simulations. 



Bo 


va 


V 


A:-grid 


V 






[G] 


[ cm • s"'] 


[cm^'' - s-'] 




[ cm - s"'] 


(kl = 2;r - 8) 


= In ■ 24) 


0.174 


1.2-10** 


1 


128' 


1.21 - 10'" 


0.86 


0.29 


0.174 


1.2-10'* 


1 


256' 


1.21 - 10'" 


0.86 


0.29 


1.74 


1.2-10" 


10 


128' 


2.9 ■ 10'" 


1 


0.37 



Notes. The outer length scale was set to Lscaic = 3.4 - 10*'cm for each simulation. The number density iij = lO^^cm"' connects the background field 
Bq with the Alfven speed = 60/ yj^mniij. All of the simulations were performed with 10^ test particles with the speed v. The n = 1 resonant fin 
values for the peaks in each simulation are shown in the two columns on the right. 



A constant absolute value of the momentum was chosen so 
that particles with a certain parallel velocity component fulfil the 
resonance condition Eq. (22 1. Consequently, a resonant value of 
J" 



to 



oj — nD. oj — nQ. 



^scale ^11 ^ 



(30) 



must be within the interval [-1; 1] for a given particle speed v. 
Since Eq. (30 1 is dominated by Q and hence Bq for particles 



propagating significantly faster than the Alfven speed, the par- 
ticle speed V has to be different for By and B^. The parameter 
setup and the resonant /i^ for each simulation is summarised in 
Table [U 

To investigate the particle scattering in different stages of 
the turbulence evolution, multiple simulations were performed. 
The most promising match between QLT and the simulation re- 
sults is expected to be in stages with small 6B/B0. According 
to this, the particle simulations were performed in the driving 
phase of the peaks as well as in the decaying phase. A simu- 
lation at maximum driven peaks would simply lead to random 
scattering where no reasonable prediction can be made by QLT. 

The simulations with the low magnetic background field Bg 
were also performed with a higher spatial resolution of 5 12^* grid 
cells, e.g. 256^ wave modes. For reasons of clarity the results of 
these simulations are presented and discussed in the Appendix 

m 



4. Results and discussion 

4.1. A toy model for wave-particle resonance 

In order to motivate our particle-scattering simulations, we want 
to present a simple model. Because of the limited applicability 
of QLT, it is necessary to understand its range of validity. In de- 
tail, the 6B/B0 ratio and the time development are of utmost im- 
portance. Different simplified scenarios are therefore presented 
below. 

We consider a circularly polarised Alfven wave with 
k - (0, 0, k\{)^ . Its magnetic field is described by 

B - Bq + 5B {e^ cos(±fe + i/') + e,. sin(±fe + 1^)) , (31) 

in the wave frame where Bq » 6B is assumed. The interaction 
between particles and the wave will then change the particles' 
pitch-angles via scatterings. A detailed description of this pro- 
cess as well as a derivation for the time dependency of jd is given 
in Appendix [a] Under the assumption given by Eq. ( 31 1. the ex- 
pression of the time evolution of fi in Eq. (|A.8[) will simplify 



A//^(f,^) = Q ^/(W^ 



6B 



cos if/ - COS [{+kvii - f2)Af + if/] 



(32) 



+kvfji - Q 

when the calculation is performed along the unperturbed orbit, 
// = constant. The change in ju is then maximal at the phase 



M 



arctan 



sin(+fcv/i - Q)Af 



1 - cos(±^VjU - 0)Af 



(33) 



and its identical solutions at periods + nn. 

To study this model, a single linearly polarised Alfven wave 
was simulated, propagating undisturbed with va - 10^ cm ■ s"' 
towards positive z-direction iw^) with a purely parallel wave 
vector k' - 2n ■ (0,0, 1)^. It should be noted that the magnetic 
field given by Eq. 



(31 



does not describe a propagating wave 
but a static disturbance. Linear polarisation was chosen to cover 
both resonances at i/i^ (respectively both signs in Eq. ( 32 1), be- 



cause interactions with /i < are caused by left-handed circular 
polarised waves and /i > by right-handed ones. The wave am- 
plitude was set to 6BIBq = 0.1, where Bq = 4.34 - 10"'* G. The 
outer length scale of the simulation cube was set to 10^ cm. 




Fig. 1. Result of a toy model wave-particle resonance after two 
gyration periods with a wave amplitude of SB/Bq - 0.1. The 
pitch-angle scattering is presented in terms of change in Afi de- 
pendent on the particles intial jUq. Each dot represents an individ- 
ual particle. The theoretical, quasilinear prediction of the maxi- 



mal change in fi by Eq. ( 32 1 deviates clearly from the observed 



scattering. This type of plot is hereafter denoted as scatter plot. 



In Gismo-Particles, 5 
position X were initialised 



10^ protons with random //, (p and 
The absolute value of v was set to 
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2 ■ 10 cm ■ s where particles for n 



-1 and n - +1 are 



resonant at jiRi - -0.28 and //r2 = +0.38 according to Eq. (22|. 

In Fig. [T[ we show the pitch-angle change of the simu ated 
particles and the QLT-predicted change as a function of the ini- 
tial pitch-angle for wave amplitude SBjB = 0.1 after two gyra- 
tion periods. The simulated particles are represented as a scatter 
plot, whereas the QLT result for different wave modes are repre- 
sented by curves. As expected, two maxima develop at the pre- 
dicted positions jiRi and ^r2- The adjoining maxima are caused 
by particles interacting with the wave nonresonantly, henceforth 
named nonresonant ballistic interactions. These decrease their 
amplitudes with increasing time. Also, their positions in fi do 
not remain constant, like the resonant interactions, because the 
balUstic maxima become more numerous and move closer to the 
true resonance. The tilt of the maxima deviates from the pre- 
diction of Eq. ( 32 1, as clearly shown by comparison with the 
theoretical curves. Furthermore, the resonances of Eq. (32i are 
located symmetrically at yU,. = ±0.33. This shift can be explained 
by the choice of the inertial system. The derived equation is 
within the wave frame. The transformation into the laboratory 
frame (under the assumption of Galilean transformation) via 



(34) 



will fit to the observed position in yu. This shift is caused by the 
magnetostatic approximation, since the difference between Eq. 
30 and its magnetostatic version (w - 0) is exactly the term 
va/v. If this ratio is small, as in most of our simulations, the 
effect is negligible. However, for the simulations with stronger 
this poses a problem, which will be discussed in the results 
section. 



The tilt cannot be explained by Eq. ( 32 1 due to its depen- 
dence on the QLT. It stems from the distortion of the orbit by 
the wave's 6B and consequently, arises by the presentation of 
Aju as a function of the initial pitch-angle, /uq. Because of the 
finite 6B/ Bo-ratio, the assumption of unperturbed orbits, which 
is the basic approximation in QLT, does not hold anymore. The 
change of fi should therefore be calculated as an average over 
the whole trajectory of each particle, which would introduce a 
nonlinear correction to the theory. A simple solution to correct 
the tilt is the presentation of Afi dependent not on the intial /io, 
but on the mean value {fj.) of the initial and the final value of the 
pitch-angle cosine. We note that the angle of the tilt does not de- 
pend on 6B/B0, but the scattering amplitude does. Consequently, 
the tilt is more visible and influences further calculations more 
as the 6B/ Bo-Tatio increases. This visibility is called effective tilt 
hereafter 

In order to prove the explanations given above, we performed 
a similar simulation with smaller wave amplitude and longer 
time development. As shown in Fig. |2] the QLT is valid within 
this scenario. As the scattering amplitude is significantly smaller, 
the effect of the tilt resulting from using the initial pitch-angle is 
negligible and the time evolution leads to a clear resonance as the 
side maxima become smaller But even under these convenient 
conditions the deviation from the QLT is still visible. 

The last case of our toy model uses an oblique wave with 
k' - In ■ (0, 1, 1)''^, which was simulated under the same con- 
ditions as before. The result in Fig. [3] clearly shows the higher 
orders of resonance with n - [-3, -2, -1, 1,2]. A prediction of 



Eq. ( 32 1 is not applied, since it would not describe the results of 



the simulation very well for two reasons. First, the assumption 
that the prediction by QLT becomes better with f ^ 00 would 
not hold because of the finite 6B/B0 and hence perturbed orbits. 
Second, the nonvanishing Bessel functions have to be taken into 



1AH+(t,<|>+) 
A[i'{t,»|/"+7r) 




Fig. 2. Result of a toy model wave-particle resonance after ten 
gyration periods with a wave amplitude of SB/Bq = 0.001 . Each 
dot represents an individual particle scattered from its initial yuo 
by AyU. The prediction by Eq. ( [32] i has been transformed into the 
lab frame and now describes the observed scattering very well. 




Fig. 3. Result of a toy model wave-particle resonance after 50 
gyration periods with an oblique wave with an amplitude of 
SBjBo - 0.001. As expected, also higher orders of resonances 
with n - [-3, -2, -1, 1,2] also developed. At this timestep Eq. 



( 32 1 would again deviate significantly, because only few parti- 
cles would fulfil the resonance condition, where others are still 
disturbed by the finite distortion of Bq and thus interact nonres- 
onantly, which produces the finite ballistic background. 



account and would modify Eq. ( 32 1. But although the ballistic 
maxima become significant in number and generate a nonreso- 
nant background scattering, the resonances are clearly visible. 
Another interesting fact is that the scattering caused by reso- 
nance at \n\ - 2 exceeds the ones at |n| - 1. This means that in 
the case of an oblique wave, higher order resonances can become 
even more important than the fundamental ones. The time devel- 
opment for the purely parallel wave mode at 50 gyration periods 
showed expectedly only the \n\ - 1 resonances (not shown in 
Fig.[3]for reasons of clarity). 



4.2. Results of particle scattering in turbulence with amplified 
modes 



We first present selected results from Lange & Spanier ( 2012|l 



which are the base of our particle simulations. However, we do 
not discuss the evolution of the peaks in detail, but focus on the 
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Fig. 4. Two-dimensional magnetic energy spectra of both peaks 
in the simulation with Z?g = 0.174 G. These are the base sce- 
narios for the particle simulations during the decay stage of the 
peaks. Therefore the left figure shows the state for the fc^ -In-Bi 
peak at f = 17 s. The right figure shows A;^ = 27r ■ 24 at f = 5.1 s. 
The colours indicate the logarithm of the total spectral energy. 



stages of the driving and decay of the peaks on which our particle 
simulations were performed. 





20 



Fig. 5. Two-dimensional magnetic energy spectra of both peaks 
in the simulation with Bq = 1 -74 G. The decay stage of the peaks 
is presented. The left figure shows the state for the peak at k'^^ - 
27r ■ 8 at f = 13.6 s, whereas the right one shows k'^^ - 2jt ■ 24 at 
t = 2.04 s. The colours indicate the logarithm of the total spectral 
energy. 



In Figs.|4]and|5]we show the situation for the decaying peaks. 
The timesteps of these figures are the starting points of our par- 
ticle simulations for the decay stage. We note that the timesteps 
were chosen to retrieve roughly the same order of magnitude 
in 5B/Bi) at the peak positions, which are 10"^. All of the peaks 
show a strong perpendicular evolution at this stage. Furthermore, 
higher harmonics of k'^^ = 27r ■ 8 at positions 16, 24, 32, 40 are 
present. These harmonics, however, do not interact dominantly 
with the particles as they have already lost most of their energy 
compared to the maximum driven stage (not shown here,) as in- 
dicated by the logarithmic colour scaling. 

As shown in the previous section on the toy model, the res- 
onance takes several gyration periods to develop. On the other 
hand, if the particles are simulated for too long, the interactions 
with the turbulence will smooth the resonances and lead to un- 
structured scattering. This is again caused by the finite perturba- 
tion of the particle orbits, which are assumed to be negligible for 
the applicability of QLT The resonances would not converge to 
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Fig. 6. Time evolution of the pitch-angle scattering coefficient 
Daa- The setup used is the low Bq, 256^ gridsize, peak posi- 
tion k!. - 2n ■ 24, driven stage. The vertical lines represent the 



predicted positions of the resonances according to Eq. ( 30 1 



a 5-function as predicted by QLT, but be scattered ballistically 
at the finite distortions of the fields. We find that a reasonable 
timescale is between 10 and 30 gyration periods for parame- 
ters used in this paper For weaker turbulence, even results at 
50 gyrations might yield a structure, but in most cases the res- 
onances develop clearly in shorter timescales. The pitch-angle 
scattering coefficient D„„ was calculated for different stages of 
the evolution. In Fig.[6]we show D„„ for the first simulation setup 
with Z?Q = 0.174 G within the 256^ grid, with the peak position 
fcy = 2;7r ■ 24 at the driven stage. We observe a convergence be- 
tween 5 and 10 gyration periods. As presented in Table [T[ the 
resonant interaction with the k'^^ = 27r ■ 24 mode is located at 
fiR = 0.29, whereas Fig. |6] shows two maxima at = 0.2 and 
fi = 0.4, which seem to move away from each other during the 
time development. 
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Fig. 7. Scatter plot for the low By, 256^ gridsize, peak position 
k'^^ - 2n ■ 24, driven stage, f = 10 gyration periods. Each dot 
represents the total change of ji of an individual particle. Three 
resonance patterns are visible and centred at fiR - -0.27, 0, 0.29. 



The real resonant structure is revealed by the scatter plot in 
Fig.|7] Indeed, the resonance is centred at the predicted position, 
but tilted for the same reasons presented with the toy model. 
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The tilt spreads the particles to a wide A/i range, which results 
in splitting of the maximum of the pitch-angle diffusion coeffi- 
cient into two maxima at both sides of the resonant pitch-angle 
— 0.29. This is because the calculation of D„„ in the QLT 



with Eq. ( 25 1 is not dependent on the sign of A//. Consequently, 



the scattering coefficient is mapped due to the square value of 
Aju to two different maxima. This also explains the movement of 
the maxima between five and ten gyration periods in Fig. |6] be- 
cause A/j increases with time. The smaller resonance at n = -1, 
i.e. jUr = -0.27, is caused by the different polarisations of the 
peaked mode. That means resonances with /i < are caused 
by left-handed circular polarised parts of the peaked mode and 
> by right-handed ones. Furthermore, in the scatter plot the 
Cherenkov resonance n = is visible. It is hardly observable in 
the Daa plot due to the dominant \n\ - 1 resonances and their 
tilt. 
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Fig. 8. Time evolution of the pitch-angle scattering coefficient 
Daa- The setup used is the low Bg, 256^ gridsize, peak position 
A:y = 2jt ■ 24, decay stage. The structure became more complex 
due to higher order resonances. The vertical lines represent the 
predicted positions of the resonances according to Eq. pO]). 



The test particle simulation with the decay phase of the 
peaked mode at A:^ = 2;7r ■ 24 shows resonant interactions beyond 
the fundamental resonance. For example, the maximum located 
near fi = 0.6 in Fig. [8] represents the n - 2 interaction. We note 
that the tilt in the corresponding scatter plot in Fig. |9] causes 
again the split of the maximum in Daa- 

As already shown by the toy model, these higher orders are 
generated by oblique Alfven waves. Those modes clearly exist 
in the decay stage of the peaked mode. For example, in Fig. |4] 
the maximum of the turbulence energy spectrum is shifted from 
the axis towards higher perpendicular wave modes. However, 
within the driven stage of the peak the energy is injected in 
purely parallel modes, which causes the dominant \n\ = 1 res- 
onances. Furthermore, the resonance pattern in Fig. |9] indicates 
that the left-handed circular polarised parts of the peaked mode 
decayed faster, which lowers the scattering frequency for par- 
ticles with fi < Q. While this observation needs further inves- 
tigation, it seems to be connected to the turbulence evolution, 
since the wave-particle toymodel does not show this behaviour 
Additionally, the magnetic background field plays a role, which 
is discussed for the results of the other parameter setup below. 
Again the Cherenkov resonance « = is visible. 
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Fig. 9. Scatter plot for the low By, 256^ gridsize, peak position 
k'^^ - In - 24, decay stage, f = 10 gyration periods. Several reso- 
nances occurred during the decay stage of the peak. The scatter 
plot reveals the resonances more accurately than the Daa shown 
in Fig.|8] especially smaller substructures are visible. 




Fig. 10. Time evolution of the pitch-angle scattering coefficient 
Daa- The setup used is the low B^, 256^ gridsize, peak posi- 
tion fcy = 2;7r ■ 8, driven stage. The resonances occurred at their 
predicted position. 



The results for the particle scattering at the k'^^ = 2;7r ■ 8 

During the driven stage 



10 



driven peak are presented in Fig. 
of the peak, the test particles interacted resonantly at \n\ = 1 
and the predicted pitch-angle jjR - 0.86. The effect of the tilt 
on the resonances is smaller compared to the k'^^ - 2n - 2A peak. 
Consequently, the maxima in the scattering coefficient are not 
split. 

The situation is different for the decay stage. The energy 
has spread significantly to oblique modes (see Fig. [4] left-hand 
frame). This leads to stronger scattering and hence a stronger tilt, 
which is clearly visible in Fig. 1 1 Even the sharp lines of the 



maximum change of /i are observed, because several particles 
were strongly scattered. In this state, Daa has no clear structure. 

The last set of simulations used an increased magnetic back- 
ground field Bo ■ Thus, the 6B / Bq ratio is decreased by an order of 
magnitude, which is more consistent with the assumptions of the 
QLT. This claim is supported by Figs. 12 and 13 where each res- 
onance has a reduced effective tilt because of the smaller scatter- 
ing frequency and represents a dominant structure. We note that 
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Fig. 11. Scatter plot for the low B'^, 256^ gridsize, peak position 
- 27r - 8, decay stage, f = 10 gyration periods. The strong scat- 
tering led to strongly tilted structures. Several particles reach the 
maximum change of /u as indicated by the sharp straight lines. 
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Fig. 12. Scatter plots for the high B'^, 256^ gridsize, peak posi- 
tion ^j'l =271-8, decay stage, f = 50 gyration periods. As in the 
toy model, the reduced SB/Bq ratio leads to decreased effective 
tilts of the resonances. After a long simulation time of 50 gyra- 
tions, the resonant interactions become very narrow. The small 
increase of the scattering near yUo - 0.5 might be caused by the 
k!. - 2n ■ 16 higher harmonic. 



the evolution time is roughly twice as long as within the weaker 
magnetic background field. A strong Cherenkov resonance is ob- 
served for the = 27r ■ 8 peak. Also both n = 1 1 1 resonances are 
visible. The small increase of the scattering rate at 0.53 could be 
generated by a resonant interaction with the higher harmonic at 
= 27r ■ 16, which is also a dominant wave mode in this stage 
(see Fig. [5] left-hand frame). 

For the k'^^ = 27r ■ 24 peak, the fundamental resonance n = 1 is 
strongest. A higher order of resonance is just slightly visible at 
fiR - 0.70. An indication of the decay of the left-handed modes 
is the absence of resonant interactions n < 0. These are visible 
during the driving of the peak, but vanish in the decay stage. 

4.3. Comparison between SQLT and particle simulations 

In the last section we present results of the SQLT approach and 
compare them to the particle simulations. For this purpose, the 
pitch-angle scattering coefficient was calculated by using Eqs. 
( 26 1 and ( 27 i for each spectrum. In contrast to the particle simu- 
lations, the calculations via the SQLT approach produce a clear 



Fig. 13. Scatter plots for the high B^, 256^ gridsize, peak A:^ = 
2;r-24, decay stage, f = 30 gyration periods. The n = 1 resonance 
is dominant, since the energy of the peaked wave mode remained 
at the parallel axis. Nevertheless, a small structure at jUq - 0.7 
indicates the n - 2 resonance. The interactions with the left- 
handed circular polarised modes vanished after the driven stage 
and are not visible during the decay anymore. 



structure in the scattering coefficient, see Fig. 14 That means 
that the resonances are not broadened by finite time and the ef- 
fect of the tilt (see the toy model section) does not influence the 
values of D„„. 
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Fig. 14. Results for the magnetostatic SQLT approach for the 
fcy = 27r ■ 24 peak within the decay stage. The calculations were 
performed for each wave species separately. The slab approxi- 
mation is comparable to the shear and pseudo Alfven waves, but 
cannot reproduce the resonance aiji - 0.58. 



As presented in Fig. [14] the magnetostatic calculations of the 
SQLT reproduce the resonant structures for n - 1,2 very well. 
However, the resonance gap for small pitch-angle cosine values 
causes a strong drop-off below < 0.25. Furthermore, a slab 
approximation was made by the projection of the total wave en- 
ergy to the parallel modes, i.e. each wave mode is assumed to be 
parallel. The slab comparison gives a good approximation of the 
scattering coefficient, except for the resonant structures caused 
by the oblique waves (|n| > 1). The scattering caused by the 
pseudo waves is comparable to the shear waves at the positions 
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of the peaks. For values of /i = [0.35; 0.55], the pseudo mode 
has higher scatter rates than the shear mode, whose coefficient 
Daa is bigger elsewhere. In the following, we restrict the presen- 
tation to the Daa sum of both wave modes, e.g. the fourth curve 



in Fig 14 
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Fig. 15. Comparison between SQLT approach and particle sim- 
ulation of the scattering coefficients Daa for the low B^, 256^ 
gridsize, peak A;^ = In ■ 8. For fi > 0.6, the SQLT coefficient is 
comparable to the particle simulations for both the background 
turbulence and the driven stage. Due to the tilt (see Fig. 1 1 1 the 
particles' Daa of the decay stage does not show a structure at the 
resonant fio. Thus, the particle simulation curve does not match 
the SQLT model nicely. 



As expected, the comparison between SQLT and particle 
simulation is best for the pure turbulent background case, as 
long as the Cherenkov resonance is small or for higher values 
of lu. Because of the resonance gap, it cannot be reproduced by 
the SQLT approach. However, the particle simulation presents 



a clear and broad interaction for n = 0, as shown in Fig. 15 
Especially for the background and driven stage, the curves in- 
crease towards jj - Q, which is caused by the Cherenkov res- 
onance. For a clearer interpretation of Daa it is again helpful 
to compare the Daa curve with the corresponding scatter plot. 
According to the previous section, the Daa is not only increased 
at yu = but also broadened by the tilt of the resonant peak. 
Nevertheless, the SQLT curve becomes comparable to the parti- 
cle simulation for ju > 0.6. Because the scattering is concentrated 
at the purely parallel peak during the driven stage, the scatter- 
ing coefficient of the particle simulation shows the resonance 
very well and is also comparable to the SQLT calculation. Even 
more complicated is the situation for the scattering at the peaked 
modes in the decay stage. The SQLT curve in Fig. 15] shows 
the resonance at the predicted position fj.^ = 0.86, whereas the 
Daa of the particle simulation seems to have no structure. This 
is again caused by the visible effect of the tilt at the resonances 



(see Fig. 1 1 



Because the effective tilt within the - In -24- peak simula- 
tions is stronger and the resonance pattern for the decay stage is 
more complex, the particle simulation results for Daa differ from 
the SQLT calculations. Nevertheless, the resonances in Fig. [T6| 
are located at the predicted positions. 

The results of the simulation with a higher magnetic back- 
ground field are presented in Figs. [17] and [18] The qualitative 
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Fig. 16. Comparison between SQLT approach and particle sim- 
ulation of the scattering coefficients Daa for the low B^, 256^ 
gridsize, peak position k'.. = 2n ■ 24. The background turbulence 
is the same as in Fig. [15] but at f = 20 gyr. The structure of the 
driven and decay stage of the particle Daa becomes clearer by 
comparison with Figs. [7] and [9] 
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Fig. 17. Comparison between SQLT approach and particle sim- 
ulation of the scattering coefficients Daa for the high B'J, 256^ 
gridsize, peak position A:^ = 27r-8. The smaller 6B/B() ratio fulfils 
the QLT approximation nicely. The particle Daa shows multiple 
maxima, where some of them are not necessarily resonances. 



comparison to the SQLT is much better with the smaller SB/Bq 
ratio. The particle scattering coefficients show clear resonant 
maxima at the predicted positions. A small shift of the SQLT 
resonances stems from the magnetostatic assumption, which sets 
cl) in Eq. ( 30 1 to zero. This is valid for w <s; f2, which is not ful- 
filled for this simulation, since = 177 s"', 1^24 - 532 s"' 
and Q. = 4264s Thus the resonances are slightly shifted to- 
wards smaller /ur, e.g. the yu^ - 0.37 for the k'^^ - 2ji -24 peak 
is yUR = 0.33 in the magnetostatic limit. Another problem is the 
resonance gap, which is significantly broader Consequently, the 
SQLT curves cannot be compared quantitatively to the particle 
simulations. 
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Fig. 18. Comparison between SQLT approach and paiticle sim- 
ulation of the scattering coefficients Daa for the high Bq, 256^ 
gridsize, peak A;^ - 2n ■ 24. Both approaches show similar struc- 
ture for the peaked modes. The shift of the resonance in SQLT 
can be explained by the magnetostatic assumption. The back- 
ground turbulence does not match well because of the large res- 
onance gap. 



5. Conclusions 

To summarise, we firstly presented a simple wave-particle in- 
teraction model to achieve a fundamental understanding of the 
underlying processes and interpretation of numerical simulation 
results. The effects of tilted and shifted resonances were ob- 
served and further compared to QLT predictions. Additionally, 
the difference between ballistic and resonant interactions could 
be shown in the context of time development. 

Afterwards, we presented three different scenarios: 

1) The first scenario used parameters to resemble conditions 
within a heliospheric plasma at three solar radii. 

2) To study the influence of the resolution, the first simulation 
setup was used with a finer grid of 512^ cells (Appendix [d|. 

3) The connection between QLT to a smaller SB/Bo ratio was 
investigated by an approach using an artificially high Bq. 

Each scenario was simulated with two different amplified wave 
modes with peak positions at A;^ = 27!- ■ 8 and 24. Furthermore, 
two different evolution stages of the peak within the turbulence 
were investigated, which allows interactions with purely parallel 
peaked modes to be studied as well as oblique modes because of 
the peak development towards a perpendicular direction. All of 
the simulations were performed with 10^ protons. Their initial 
speed was chosen to be in resonance with the peaks. In addition, 
every turbulence spectrum was used for a magnetostatic QLT 
approach to calculate the pitch-angle scattering coefficient. We 
presented resonance patterns with different methods: the scatter 
plot Afi (yU()) and the scattering coefficient (jko) for both the 
numerical and semi-analytical quasilinear approach. 

Our results show a good agreement between hybrid MHD 
particle simulations and QLT calculations for the scattering co- 
efficient Daa for a turbulent broad-band spectrum. The compar- 
ison reveals that the SQLT calculations based on the MHD spec- 
trum are missing the Cherenkov resonance (n - 0), which is 
not covered by the QLT calculations used because quasilinear 
theory produces a singularity at yu = 0. Furthermore, the reso- 



nance gap at small values of is caused by the power spec- 
trum, which does not have sufficient energy at high wave num- 
bers within the dissipation range of the turbulence. The direct 
particle simulations, on the other hand, show problems linked 
to broadened resonances: Because of the finite simulation time, 
the wave-particle resonance peaks are broadened, which in turn 
yields a significant scattering background. With increasing sim- 
ulation time, especially in simulations with higher SBjBo, the 
resonances did not become narrow but scattering would start 
to randomise. Consequently, the assumption of 5-shaped reso- 
nances does not apply to real particle scattering. Furthermore, 
the tilt in A/i, as shown in the scatter plots (e.g. Fig. |7]l, causes 
spreading of the resonant peaks over an interval of for the 
particle simulations. 

Due to the tilt and the broadening of the resonances, the 



application of the scattering coefficient resulting from Eq. ( 25 1 
is questionable. The Daa results are either difficult to interpret 
or without any reasonable structure. Consequently, the QLT ap- 
proach does not compare well to the particle simulation results. 
We point out that in these cases a very nice tool in numerical 
simulations is the scatter plot. By evaluating the total change of 
ji versus its initial state /^o, even small resonance patterns be- 
come visible. This can be used for correct interpretation of the 
scattering coefficients. Especially in cases with multiple over- 
lapping resonances (e.g. Fig.|9|, the scatter plot yields important 
information. It should be noted, however, that the analysis of 
particle-scattering data requires more in-depth research. 

A nice comparison between QLT and particle simulations 
is achieved with smaller 6B/B0 ratio of the turbulence and the 
peaked modes. This is not unexpected, because the assumption 
of unperturbed orbits is more reasonable. In this case the effec- 
tive tilt in the scatter plots decreases significantly and Daa be- 
comes more structured. Unfortunately, the resonance gap in the 
QLT approach expands in these scenarios, which leads to smaller 
absolute values of the scattering coefficient. This large gap stems 
from the used spectrum, where the magnetic background field 
was increased, leading to less resonant waves. 

We conclude that for realistic particle transport, especially 
around /i ~ 0, the QLT approach does not yield physical scatter- 
ing coefficients. This is not unexpected, as it has been discussed 
in literature akeady ( Shalchi et al.||2004) l We are not aware of 
any self-consistent theory describing the numerical particle re- 
sults derived from our hybrid simulations. Further investigations 
are needed to disentangle the different transport processes in- 
volved in order to develop a new kind of transport theory. We 
also note that even our numerical approach has limited validity 
for high SBjBo ratios. We will discuss more sophisticated nu- 
merical methods for transport analysis in a forthcoming paper. 
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Appendix A: Time dependency of the pitch-angle 
cosine 

The interaction of a charged particle with the fluctuation 6B of 
an Alfven wave forces the particle to leave its gyro orbit. This 
procedure is referred to pitch-angle scattering. Consequently the 
change of parallel momentum is connected to this process via 
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~ W- Then the time derivative of is given by the Lorentz 
force 



dt 



ec 

2En 



(A.l) 



where Ep is the energy of the particle, which is assumed to be 
constant as the scattering is purely elastic. Since the parallel mo- 
mentum remains constant during an unperturbed gyration, only 
the perpendicular momentum changes due to the circular move- 
ment - + ipy = Pj_() exp[+/(0o - i^t)]. This leads to a total 
change of fj. for a particle with intial phase (po to the magnetic 
field with 



ec 



dt 2pEp 



J](+OPxo exp[+/(0o - nt)]6Bi/r), (A.2) 



([Lee & Lerchel |1974|l A similar approach is given by 
Schlickeiser (2002 ) using Eq. 19 In this case the solution along 
the characteristics of the generaUsed force term - fi\s given 
by Eq. (12.2.4b) Schlickeiser '{ 2002") l This result reduces to Eq. 
A.2 by using the magnetostatic approximation, which assumes 
the electric field fluctuations to be negligible, SE - 0. The co- 



ordinates used are still Eqs. 20 The fluctuation 5B, i.e. of an 
Alfven wave as used in the presented case, is given in Fourier 
space by 



5B,,(r) = -^ f d"A:(5B,,(A:)exp(/A:r). (A.3) 
on^ J 

The exponential function can be described by the generating 
Bessel functions /„ 

DO 

exp (ikr) = Jnix) exp /(fcy vyf + - 0,) + Or)), (A.4) 

«=-oo 

where the argument of J„ is 

y a/1 — i]^ 

(A.5) 



Q 



and *P = cot ^{ky^lk^), kj_ - + k^.. The total time derivative 

of the pitch-angle cosine, separated in parallel and perpendicular 
interactions, reads then 



8m 



d^i 



■ ^(+0/'±o exp[-i-i0o] 



dt 2pEp ^ 

J d^k 5B,/k) ^ y„ (z) : 



exp [in{^ - 0o) + it{k\\V\\ ± O - nQ.)], (A.6) 

For Alfven waves 5B is aligned towards k x e,BQ and conse- 
quently 6Bi f {k) = 6B(k){+i) exp(+/*l'). The time integration and 
the identity 

—J„iz) ^ Jn^iiz) + J„-dz) (A.7) 



( Abramowitz & Stegun 1965|l| then gives the time dependency of 
the pitch-angle cosine with 



Ap(t) = fi(t) - fio 



ecpj_o 



d^k 6B(k) X 



2pEp 2,7, 

Un+l(z) + Jn~\(Z)\ X 

/I--CXD 

exp [;■«(*!'- (^o)] f df' exp [;f'(A:||V||o - nO)]. (A.8) 
Jo 



The application of a purely parallel propagating wave, as shown 
in our toy model, will then simplify this equation because 
z{ki_ = 0) = 0. Thus, all Bessel functions vanish, except 7o(0) = 
1. This is the case for n - +1. Furthermore, the k integration re- 
duces by the assumption of a single wave k - 5(k\\ - ko)S(kj_)k, 
which leads to the presented form in section [4T| 



Appendix B: Derivation of the pitch-angle diffusion 
coefficients 

In the simulations, the turbulence consists of Alfven and pseudo 
Alfven waves, thus the pitch-angle diffusion coefficient must be 
calculated separately for the two modes. The two modes are 



decomposed using the method presented Maron & Goldreich 
(2001) For the pitch-angle diffusion coefficient for Alfven 
waves, |Schlickeiser ( 2002j i gives 



2Q2(1 - ^P-) 



B^ 



2 



kHiXcS) 



1 _ 

k\\v\ 



nJ„(vj_kj_/Q.) 



v±kj_/Q. 



P„,A(k), (B.l) 



where p, v and Q again are the particle's pitch-angle co- 
sine, speed, and gyrofrequency, Vj_ - v -y/l -p^, Pxx,a the xx- 
component of the Alfven mode turbulence power spectrum ten- 
sor, and Bo the background magnetic field. For the resonance 
function K (k, oj) we assume no damping, which gives the delta 
function '??(k, w) - jT6{k\\v\\ - a> + nO.). Thus, in the magneto- 
static limit {(JL> - 0) we have 



^ I k TiS {ki\vii + riD.) X 

1— 



nJn(v±k±/Q.) 



PxxA^). (B.2) 



For the magnetosonic waves, [Schhckeiser ( 2002|l| gives for the 
fast mode on the cold plasma limit 



D 



2Q?(l -p^) 



OO ^ 



k'R{k,Lo)x 



(B.3) 



for high particle velocities, <s; v. As the polarisation of the 
fast mode on the cold plasma limit is the same as for the pseudo 
Alfven wave, it is straightforward to show that the same form 
also holds for the pseudo Alfven waves, with the difference in the 
dispersion relation. However, by using the magnetostatic limit 
with (L) = and again assuming no damping, we get for the 
resonance function Kik, a>) = 7r(J(fe||V|| + Q.) , thus giving for the 
diffusion coefficient 



2^2(1 _^2) 



B^ 



oo ^ 



k 7t6 {k\\v\\ + nQ) X 



(B.4) 



Unlike for the Alfven waves, the Cherenkov resonance, 
« = 0, is nonzero for the pseudo Alfven waves and has to be 
considered separately. In this case, the resonance function is 



f^n=0 - 



;r5(v||) 



(B.5) 
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This results in 



For the Alfven waves, using the same discretisation, we get 



L>..„ „ - ^ (V||) X 



J flf3A:^f/^(v^^^/Q)]'p„,Kk), (B.6) 

which has a singularity at vy - 0, and equals zero elsewhere. 
Consequently, this term is not used by our model. 

For other terms in the sum over n, we again use the magneto- 
static approximation u -Q, thus giving the resonance condition 



Hik, u) - nS {k\\v\\ + nQ) 
and the pitch-angle diffusion coefficient 

2Q2(1-|/2)^ f 



(B.7) 



•>l 



2/ 



d kn6{k\\v\\ + nQ) x 



(B.8) 



Appendix C: Discretisation of the pitch-angle 
diffusion coefficients 



For the numerical calculation of the pitch-angle diffusion co- 
efficient, we consider the spectrum to be continuous in par- 
allel direction, but discrete in the perpendicular direction, i.e. 
Pxxikxky, k,) - P„(hAk, iAk, L), with h, i - -M, . . . , M. In this 
manner, the integrals over k^ and ky can be written as a sum, 
while the integral over k, is evaluated using the delta function. 
Then, the diffusion coefficient for the Pseudo- Alfven waves can 
be written as 



,n*0 



Blh\ 

2 2 Ak^[J'„(v^k^/n)fP,,AhAk,iAk,nil/v^^). 

n=-oo hJ=-M 

(C.l) 

In this equation, (nQ.)/(p.v) represents the parallel wavenum- 
ber at which <5(A:||V|| + «Q) is nonzero. Consequently, as turbu- 
lence data is available only at discrete wavenumbers, we define 
(nf2)/(juv) - IAk. Thus, with I restricted to integer values, we 
find the values of (vy - fj.v) at which the value of can be 
solved. For a particle with v = SljimAk), with m < I < M, the 
pitch-angle is given as /i = m/l. Applying this discretisation, we 
have 

V|| - v(m/l), 
Vj_ — V ■sjl - nP-/P-, 
and k^ = ^kl + kj = Ak ^li^ + fi (C.2) 



and thus 



m 

M//,«#0 M 

z z 



I X 



j: 



n=-M/l h,i=-M 

Ak^ P^,p{hAk,iAk,nlAk) 



(C.3) 



/m\ / m \ 

Z),M(7) = 27rQjl--^jZx 



Mil M 

Z Z 

n=-M/l li,i=-M 



Ak^ P„A{hAk,iAk,nlAk) 



(C.4) 



Appendix D: 512^ Results 

In this section we present the results of the investigation of the 
resolution by using a spatial grid of 512^ cells. Again, the MHD 
simulations were performed, first for the background turbulence, 
afterward with the peaked modes at = 2;r ■ 8 and 24. The 
results are in conformance with the other setups. In particular, we 
observed the generation of higher harmonic wave modes again. 





20 40 60 



20 40 60 SO 



Fig.D.l. Two-dimensional magnetic energy spectra of the decay 
stage for both peaks in the simulation with - 0.174 G and 
higher resolution with a grid of 512^ cells. The left figure shows 
the state for the peak at k'^^ -2ji-%sit - 14.45 s. The right figure 
shows the k'^^ -2ji-2A peak at f = 3.4 s. The larger Fourier space 
grid reveals the higher harmonics of the k'.. = 2;r • 24 peak. 



These harmonics are also visible for the peak at k'^^ - In ■ 24 

within the higher resolved simulation with a 512^ grid. In this 
simulation setup the number of active modes is eight times 
greater, which means the antialiasing edge is shifted by a fac- 
tor two to k' - 2-n ■ 86. Consequently, the higher harmonics at 
k'^^ - 27r-48 and 72 are visible. The evolution of the peaks in both 

simulations is comparable to the 256-' grid simulations. A domi- 
nant energy transport towards high perpendicular wavenumbers 
is observed. 

The test particle simulations led to comparable resonance 
patterns. As presented in Fig. |D.2| the resonant interactions at 
//r = 0, 1 and -0.95 are strongly tilted and a significant amount 
of particles reach the maximum A/i, as indicated by the sharp 
thresholds. The scattering coefficient is in this case again with- 
out any structure and hence not shown here. The stronger scat- 
tering is primarily caused by the higher wave modes, which are 
not truncated by the antialiasing anymore and hence contribute 
to the wave-particle interactions. The test particle simulations 
were performed for the decay stage of the peaks only because of 
their huge computational effort. 
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tion = 27r ■ 8, decay stage, f = 25 gyration periods. The higher 
resolution of the grid causes slightly stronger scattering, due to 
the higher amount of active modes. The resonance patterns are 
comparable to the 256^ grid. 
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Fig.D.3. Scatter plots for the low Z?g, 512^ gridsize, peak posi- 
tion = 2:7r-24, decay stage, t -25 gyration periods. Because of 
the higher resolution, higher harmonics of the k'^^ - 2jt ■ 24 peak 
have also developed and interact with the particle. This leads to 
more diffuse scattering. The resonances are barely visible. 



Additionally, an increase of the scattering rate at the k'^^ - 

2n ■ 24 peak within the 512^ grid was observed. Consequently, 
the resonances are not as significant as in the smaller grid. When 
comparing Fig. [9|and |D.3| it is harder to recognize the resonant 
structures. 

As discussed previously, the 512^ gridsize particle simula- 
tion within the decay stage led to strong effectively tilted res- 
onances and consequently a very unstructured D„„ curve. This 



can be observed in both Figs. D.4 and D.5 Thus, the results dif- 



fer greatly from those given by SQLT. 
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